Lyapunov exponents and transport in the Zhang model of Self-Organized Criticality. 
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We discuss the role played by the Lyapunov exponents in the dynamics of Zhang's model of Self- 
Organized Criticality. We show that a large part of the spectrum (slowest modes) is associated with 
the energy transport in the lattice. In particular, we give bounds on the first negative Lyapunov 
exponent in terms of the energy flux dissipated at the boundaries per unit of time. We then establish 
an explicit formula for the transport modes that appear as diffusion modes in a landscape where the 
metric is given by the density of active sites. We use a finite size scaling ansatz for the Lyapunov 
spectrum and relate the scaling exponent to the scaling of quantities like avalanche size, duration, 
density of active sites, etc ... 

PACS number: 02.10.Jf, 02.90+p, 05.45. +b, 05.40. +j 

I. INTRODUCTION. 

Within the last 10 years the notion of Self-Organized Criticality (SOC) has become a new paradigm for the 
explanation of a huge variety of phenomena in nature and social sciences. Its origin lies in the attempt to explain 
the widespread appearance of power-law like statistics for characteristic events in a multitude of examples like the 
distribution of the size of earthquakes, 1/f-noise, amplitudes of solar flares, species extinction .... to name only a few 
cases • In this paradigm, the dynamics occur as chain reactions or avalanches. Furthermore, a stationary regime 
is reached, where the average incoming flux of external perturbations is balanced by the average outgoing flux that 
can leave the system at the boundary, or by dissipation in the bulk and there is a constant flux through the system. 
In this stationary state, referred to as the SOC state, the distribution of avalanches follows a power law, namely there 
is scale invariance reminiscent of thermodynamic systems at the critical point. A local perturbation can induce effects 
at any scale and there are long-range spatial and time correlations. In other words, in this paradigm the system 
reaches spontaneously a critical state without any fine tuning of some control parameter. 

Several models have been proposed to mimic this mechanism, like the sandpile model the abelian sandpile 
[^j or the continuous energy model 0]. The results available are mainly numerical and only a few rigorous results 
are known. The numerical simulations report the following behaviour. Fix an observable, say x, measuring some 
property of an avalanche (duration, size, etc ...), and compute the related probability Pl{x) at stationarity for a 
system of characteristic size L. The graph of Pl(x) exhibits a power law behaviour over a finite range with a cut-off, 
corresponding to finite size effects. As L increases the power law range increases, leading to the conjecture that a 
critical state is indeed achieved in the thermodynamic limit, namely that Pl{x) behaves like as L — > oc. t x is 
called the critical exponent for the observable x. There is apparently no control parameter to tune to attaining the 
critical state. Despite the large number of papers written on the subject some basic problems arc still open. 
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Guided by the wisdom coming from renormalization group analysis and phase transitions in equilibrium systems, it 
seems natural to look for a possible classification of the models into universality classes characterized by a set of critical 
exponents, for a family of "relevant" avalanche-observables. However, the link between the "criticality" of the "out 
of equilibrium" SOC models and the usual statistical mechanics of phase transitions in equilibrium systems remains 
to be clarified ||. Furthermore, apart from the fact that the commonly studied observables (size, duration, area, 
giration radius) do not necessarily consitute a complete set allowing one to classify the models, even the computation 
of the critical exponents t x from numerical data is not easy and there is no agreement yet on the way to do it. It is 
clear that the simple measurement of the slope of Pl{%) in the linear range of a log-log plot is not reliable, due to 
the finite sample fluctuations and because the explicit form of the cut-off is not known in general. The computation 
of t x from the behaviour of the moments is certainly a better way. However there is no agreement yet wether one 
should use a finite size scaling treatment or more sophisticated methods (like multifractal analysis Q). Therefore, 
at the moment, the identification of a (supposed) universality class seems problematic. Finally, the central question 
is: what exactly does the knowledge of the critical exponents t x teach us about the model.?. 

An alternative approach to better understand the behaviour of SOC models can also consist of studying the 
microscopic dynamics and to infer information about the macroscopic behaviour from this analysis. A detailed analysis 
can, at first sight, seem useless since the conventional wisdom from classical statistical mechanics is that microscopic 
"details" are irrelevant, and only structural properties like conservation laws and symmetries are essential. However, 
as mentionned above, the theory of SOC has not yet reached the level of understanding of classical critical phenomena. 
It suffers in particular from the lack of a thermodynamic formalism and notions like Gibbs measures and free energy 
can a priori not be used. On the other hand, by having a precise description of the dynamics of the finite size system 
one can expect a better understanding of the thermodynamic limit and decide which components in the models 
definition are really "relevant" , and which information docs the usually computed quantities (like critical exponents) 
actually give us. 

This is the essence of the program we developped in J|-|ll|. We found that Zhang's model of SOC H can be 
fruitfully studied with the tools of hyperbolic dynamical system theory. Then we were able to extract unexpected 
results, establishing in particular a formula relating the critical exponent of avalanche size to the spectrum of the 
Lyapunov exponents. In this paper we develop this point of view and make a further step towards unterstanding 
the dynamical properties of this model and their link to the SOC state. We first define the model as a hyperbolic 
dynamical system of skew-product type. We then define in this setting two different time scales : the local time 
which is the natural time for the dynamical system, and the avalanche time, related to the avalanche duration. We 
introduce a natural invariant measure to characterize the statistical properties at stationarity, and we relate the 
avalanche observable statistics to the ergodic local time average. We then discuss the role played by the Lyapunov 
exponents in the dynamics and their relation to the energy transport and to the average avalanche observables. We 
show that random excitation induces a positive Lyapunov exponent, while the relaxation dynamics corresponds to 
negative exponents. Furthermore, we show that a wide part of the spectrum (slowest modes) is associated with energy 
transport in the lattice. In particular, we give bounds on the first negative Lyapunov exponents in terms of the energy 
flux dissipated at the boundaries per unit of time. We establish an explicit formula for the transport modes, that 
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appear as diffusion modes in a landscape where the metric is given by the density of active sites. Except for the first 
modes, they differ dramatically from the normal diffusion modes one would obtain by assuming a uniform density of 
active sites. It has been argued in Jl2| that SOC requires a wide separation between the excitation and relaxation time 
scales (slow driving). We actually show in this paper, as a consequence of our more general analysis, that the dynamics 
of the Zhang provides naturally this separation, and that the infinitely slow driving limit is actually reached as the 
size of the system goes to infinity. We then show that, using a finite size scaling ansatz for the Lyapunov spectrum, 
one can relate the obtained scaling exponent to the scaling of quantities like avalanche size, duration, density of active 
sites, etc .... 

II. DYNAMICAL SYSTEM DEFINITION AND BASIC PROPERTIES. 

A. Definition. 

The Zhang's model B, widely inspired from the Bak-Tang-Wiesenfield precursor model 0, has been introduced as 
a possible example of a model which "self-organizes" into a critical state in the thermodynamic limit, namely without 
fine-tuning of a control parameter. Its beauty lies in its simplicity. 

Let A be a d-dimensional sub-lattice in 2Z d , taken as a square of edge length L for simplicity. Call N = #A = L d , 
and let 9 A be the boundary of A, namely the set of points in 2Z d \ A at distance I from A. Each site i € A is 
characterized by its "energy" X,-, which is a non-negative real number. Call X = {Xi} ieA a configuration of energies. 
Let E c be a real, strictly positive number, called the critical energy, and A4 = [Q,E C [ N . A configuration X is called 
"stable" iff X S M. and "unstable" otherwise. If X is stable then one chooses a site i at random with some probability 
and add to it energy <5, where 6 is set to I in this paper (excitation). If a site i is overcritical or active (Xi > E c ), 
it loses a part of its energy in equal parts to its 2d neighbours (relaxation). Namely, we fix a parameter e G [0, 1[ 
such that, after relaxation of the site i, the remaining energy of i is eXi, while the 2d neighbours receive the energy 
^-^ei ■ Note therefore that there is local conservation of energy. If several nodes are simultaneously active, the local 
distribution rules are additively superposed, i.e. the time evolution of the system is synchronous. The succession 
of updating leading an unstable configuration to a stable one is called an avalanche (a more precise definition of an 
avalanche will be given below). There is dissipation at the boundaries namely the sites of 9 A have always zero energy. 
As a result, all avalanches are finite. The addition of energy is adiabatic. When an avalanche occurs, one waits until 
it stops before adding a new energy quantum. Further excitations eventually generate a new avalanche, but, because 
of the adiabatic rule, each new avalanche starts from only one active site. Note that relaxation depends on local 
conditions while excitation is conditioned by global constraints (all sites are quiescent). It is conjectured that a critical 
state is reached, independently of E c , at least for large E c values [] 



1 Strong deviations from a power law have been observed for small E c in one dimension 
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B. The Zhang's model as a dynamical system. 



Because all avalanches are finite (for finite L), and since we are not interested in the transients, one can, without 
loss of generality take all initial energy configurations XeM. All trajectories starting from M. belong to a compact 
set B. Call M. = B \ M. . M contains the set of all unstable energy configurations achievable in an avalanche starting 
from an energy configuration in A4. 

Fix e > 0, and call a — ^r- Let h be the Heaviside function. Define H : R N — > {0, 1}^ such that i?(X) is the 
vector {h(Xi)} i=1 N . Call X c the vector {E c } i=l N . Finally, let A be the discrete Laplacian. The dynamics is 
defined by the mapping F : B — > B such that: 

F(X) = X + aA [H(X - X C ).X] (1) 

which redistributes the energy of the active sites in equal parts to the neighbours after one relaxation step. Note that 
F is the identity if no site is active, i.e. if X E M, and that it is piecewise linear, (i.e. linear on sub-domains Bk E B). 
F is a (singular) diffusion operator and a the diffusion coefficient. 

It is useful to encode the dynamics of excitation in the following way. Let be the set of right infinite sequences 
a = {a\, . . . , afe, . . .} , afe E A, and a be the left shift over E^, namely era = 0203 .... The elements of T,^ are called 
excitation sequences. The set f2 = x B is the phase space of the Zhang's model and X = (a, X) is a point in f2. 
The Zhang's model dynamics is given by a map of skew-product type F : ft — > f2 such that: 

XEA4^F(X) = (a.a,X + e a ) (2) 
X E M => F(X) = (a,F(X)) (3) 

The knowledge of an initial energy configuration X, and of an (infinite) sequence of excited sites a (resp. of an 
initial condition X) fully determines the evolution. One can give a probability distribution i>l corresponding to 
a random choice for the excited sites. In the original Zhang's model the excited sites are choosen at random and 
independently with uniform probability. This corresponds to giving the uniform Bernoulli measure. Throughout 
this paper we will often think of the left Bernoulli shift on as represented by the system z — > N.z mod 1, z E [0, 1]. 

In the following we will denote the two projections on the first and second coordinate by 7r u (X) = a, and 7r s (X) = X. 
The superscripts u, s mean respectively unstable and stable and correspond to the expansion (resp. contraction) prop- 
erties of the dynamics. Let DFj be the tangent map of F at X and DF*^ the t-th iterate. As shown below tt u (DF-^) 
is expansive whereas 7r s (_DF-£.) induces contraction. In the following we will use the notation X(i) = F'(X) (resp. 
X(i) = tt s (F'(X))). Furthermore note that ir s (DF ± ) = DF X , and that £>F X = /, the identity matrix over IR^, if 
X E M.. 



Consider a point X E f2. Its trajectory is intermittent, composed of bursts of excitation of the sites 01, 02, ■ ■ ■ a„, 
for those times t such that X(t) E M, followed by relaxation periods when X(i) E M. Define the following hierarchy 
of waiting times: 

70 (X) = (4) 
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(Ti(X)=Jnf {X(t)<EM}, »>1 (5) 
7i (X) = inf {X(t) G.M}, i>l (6) 

For i > 1, (Tj(X) (resp. 7i(X) ) is the starting time (resp. ending time) of the i-th avalanche occuring during the 
evolution of X. Therefore the avalanche duration of the i-th avalanche is : 

r i (X)= 7< (X)- - i (X) (7) 

In the same way, one defines : 

w i (X) = o- i (X)-7 i _ 1 (X) (8) 

which is the number of excitations between the end of the avalanche i — 1 and the beginning of the next avalanche. 
In this way, one introduces naturaly two time scales : the local time t corresponding to one step of iteration in the 
dynamics, and the avalanche time t% corresponding to the duration of an avalanche (a similar description has been 
used in Hi). 

The waiting times are useful for defining the usual avalanche observables. The number of relaxing sites for a given 
configuration is : 

r(X) =#{jeA, Xi> E c } (9) 

The avalanche size is 

r(X) 

,s(X)=^r(F t (X)) (10) 
t=i 

where: 

r(X) = inf{F t (X)GX}-l (11) 

is duration of the avalanche occuring when exciting the site a± in a stable energy configuration X. It is zero if one 
drops energy without relaxation. 

The structure of an avalanche can be encoded by the sequence of active sites A(X) = |yl t (X) j „ where 

A t (X.) = {j G A\Xj(t) > E c }. (Note that Ai(X) is non empty and equals {a±} iff X + e ai is active). Correspondingly, 
there exists a partition^ of x A4 into domains Vi t h — [i] x -M-i,k where [i] is the set of sequences in having a 
first digit i, such that for any energy configuration X G M.i t k the excitation of the site i leads to the same avalanche 
(namely the same sites relax at the same time). Under some moderate assumptions (see fll|| ), this allows us to define 
a symbolic coding for the avalanche and a transition graph giving the transition rules between successives avalanches, 
and to show that the dynamical system admits a unique, fractal, invariant set. The boundary of the domains 7\fc 
constitutes the singularity set of F, called S. This is the set of points where F is not continuous. 



2 This partition is induced by the partition of B into domains of continuity for F 
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C. Stationary state and probability of avalanche observables. 




Let fiL be an invariant measure for the dynamical system |fi,F|, where L refers to the lattice size, namely 
lil(F^ 1 (A)) = lil{A) where A E f2 is a measurable set. Since f2 has a product structure where the unstable foliation 
is always transverse to the stable one, and since the dynamical system is a skew product, {ml — vl x Mil where 
vl is the induced measure on the unstable direction or excitation measure, and lll is the induced measure on B or 
measure on the energy configurations. For simplicity we will assume that vl is a Bernoulli measure, namely that the 
successive excited sites are chosen independently with fixed rates. Once we have fixed the distribution of excitation, 
we are interested on the possible lil measures. Of special physical importance are the measures obtained by iterating 
the Lebesgue measure /iLe60 on A4, that is burgoo i Yn=o ^(^l x ULeb)- When the excitation measure vl is itself 
the Lebesgue measure on [0, 1] (corresponding to choosing the excited sites with uniform probability) the measure 
obtained is called the Sinai- Ruelle-Bowen measure (SRB). More generally, we will call (conditional) SRB the measure 
linin^oo i X^Lc) 1 ^( v l x LiLeb), for & fixed Vl. This is a natural invariant measure from the physicists point of view 
since it gives the ensemble average with respect to typical initial energy configurations. 

It is common in the SOC literature to assume ergodicity. In our setting, the physically relevant ergodic property 
is equivalent to assuming that the SRB measure is unique. Proving the ergodicity in the Zhang's model, is clearly a 
difficult task which is beyond the scope of this paper. We note however that this point has been widely discussed in 
a previous paper [Tl| | , where strong mathematical arguments in favour of this were given. Actually, ergodicity was 
proved, but restricted to the one dimensional model, and to some E c interval. A general proof is under construction 
and will be published elsewhere On physical grounds, note a contrario that the failure of ergodicity would lead 
to a stationnary state depending on initial conditions. This would contradict the implicit SOC assumption that the 
stationnary state is unique. In the following, we will therefore assume that ergodicity holds and that {ml is the unique 
SRB measure. This implies in particular the almost-sure equality between the ensemble average and the time average, 
namely, if <j> is some observable, (a function Q — > IR, integrable with respect to {ml) 



for a typical (namely Lebesgue almost surely) initial condition X. Here, and in the following l will denote the time 
average, while El[] will be the ensemble average, on a lattice of size L. 

From the dynamical systems point of view lil is the natural object to deal with. However, in the SOC literature 
one is more interested in the probability distribution of some avalanche observable and its scaling properties in the 
thermodynamic limit. Fix an avalanche observable, say s. Call V s the union of domains V^k such that the avalanche 
corresponding to each domain V^k has the same size s. Then the probability of having an avalanche of size s, by 
excitation of a stable configuration, is Prob[s(jC) = s|X E V] — j$pj — ^jja) ■ ^ n this definition we include the 



3 Or any absolutely continuous measure, which corresponds to selecting the initial energy configuration with a probability 
distribution having a density. 




(12) 
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avalanches of size zero (excitation without relaxation). However, it is more natural from the SOC point of view to 
exclude this case. We therefore define Pl{s) as the probability of having an avalanche of size s strictly larger than fiFl. 



p L[8] 6* MPs) 



S > 1 



PL 



(13) 



where pl = f Prob s(X) > 1,X 6 M. is the probability of initiating an avalanche. The average with respect to Pl[s], 
denoted further on by <>l, is : 

(I 



(14) 



where tp is some real function, and ££ is the maximal value that the observable s can have on a lattice of size L 
(note that depends also on E c ,e,d but is bounded if L < oo). The same definition holds for any other avalanche 
observable. From the ergodic theorem: 



1 - 

< ifr(s) >l= lim - y^ip(si) 

n — >oo n — ' 



(15) 



where Si is the size of the i-th avalanche occuring in the trajectory of a generic point X. 
One has: 



PL = Vl 



N 



(J {ai=i,X t e[E e -l, E c [} 



N 



where: 



p L (i) = v L {j).n L {X t G[E C -1, E c [} 



(16) 



(17) 



is the probability that an avalanche starts at the site i. Note that the probabilities Ph{i) depend a priori on i even if 
the excitation measure is uniform. In this case, however, ( |l6| ) reduces to 

N 



PL 



-^■J2vl {X, e [E C -1,E C [} 



(18) 



Fix X and T, then call n(T, X) the number of complete avalanches occuring until local time T for the initial 
condition X. Obviously, n(T, X) — > oo as T — > oo, VX. Then from the ergodic theorem : 

n(T, X) 



Pl = lim 



T^oo T 



(19) 



n(T,X) n(T,X) 

One can decompose T has : T = Tj + u>i + K(X.) where K(X.) is some residual time, finite, whatever 

i=l i=l 

T, whatever X (K(X.) is bounded by the largest avalanche duration). Note that r,; (resp. u>i) stands for Tj(X) (resp. 
£jj(X)) but we removed the X dependence in order to simplify the notations. Then, as T goes to infinity: 



4 In view of the expected critical behaviour as L — » oo, one usually writes a scaling form Pl(s) = ^fvy where /l(s) is a cut-off 
term accounting for finite size effects on large scales. Pl(s) is not denned for s = unless assuming very special properties for 
his). 
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n(T,X) n(T,±) n(T,X) _ n(T,X) E?iT' X) ' 

T ~ v-«(T,X) V n(T,X) v-n(T,X) ^n(T,X) T 

Call: 

*(T,X) 



= T lim E w, = Mi(A<) (20) 



def 

i=l 



the probability of dropping energy in the system, at a given time, (the equality holds for /z^ almost-every X from the 
crgodic theorem). cjf,(i) = Pro6[ai = i,X € At] = is the probability of dropping energy on the site «, at a 

given time, and is called the driving rate in the literature |l5| ]. One has: 

PL = = (21) 

< T > L < T > L 

where < r >l is the average avalanche duration. 

There exists an important relation linking the avalanche averages (average with respect to Pl) to the local time 
average (average with respect to fiL,)- Let <fi : Q — ► R be some observable such that 0(X) = whenever X 6 M. 
A related avalanche observable can be defined by summing the values that <f> takes in one avalanche. Namely, call 
/i(X) = X^ X |x) ^(-^(*))' (^ n important example is when 4>(K) — r(X), the number of active sites in one step. 
Then /i(X) is the size of the i th avalanche in the trajectory of X). One obtains: 

n(T,X) 7i (X) 

^ = ^t E E *(*(*)) 

«=1 t=<T;(X) 

which yields: 

H=PL-{f) L (22) 

In particular : 

rL=Pi.(s) L (23) 
Finaly we define the probability that a site i is active (often called the density of active sites in the literature : 

p L (i) d = » L [Xi>E c ] (24) 

and 

1 W 

pf^ is believed to act as an order parameter in the Zhang's model. 



3 We will keep this terminology throughout this paper though p_l(i) is not strictly speaking a density since X/^-i A^W ^ 1- 
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III. DYNAMICAL PROPERTIES AND LYAPUNOV EXPONENTS. 
A. Jacobian matrix and Lyapunov exponents. 

Due to the piecewise affine structure of the map F, the Jacobian matrix D Fx plays a central role in the Zhang's 
model, since it characterizes the energy transport. Indeed, the entry DF^ - is the ratio of energy flowing from site 
j to site i in t times steps for the initial condition X. Define Zk(X(t)) — H(Xk(t) — E c ). This is a random variable, 
taking value if Xk(t) is stable, and value 1 otherwise, whose probability distribution is induced (at stationarity) 
by the invariant measure More precisely, Prob{Zk(X(t)) = 1] = piik). Let Z(X) = {Zk(X)}^ =1 , and call 
S(X) = AZ(X)7 (equivalently S(X) is the matrix of entries Sij(X) = AyZj(X)). S is the "toppling" operator of 
the Zhang's model. The jacobian matrix writes DFx — I + a.S(X), while DF^ is given by : 
t 

DF^=I + aJ2S(X(t )) + a 2 ]T S(X(*i))5(X(to)) + . . . (26) 

to=l t>t 1 >t a >l 

+a r ]T 5(X(t r _!))5(X(t r _ 2 )) . . . S(X(t Q )) + ... + a*S(X(t))S(X(t - 1)) . . . S(X(1)) 

t>t r _l>t r _2--->*0>l 

Therefore, the generic term (say of order r) is a a "propagator" transmitting the energy in r times steps. Note that 
this formula is exact. It calls for the following remarks: 

• The maps 5(X) do not commute, and they depend on the state. This is a key difference from the Dhar's model 
since it induces a non abelian structure and a "toppling" operator depending not only on the site, but also on 
the whole energy configuration. In particular the propagator is not a mere polynomial of the Laplacian. 

• The evolution depends a priori on the whole trajectory and therefore the strong memory effects expected in a 



critical phenomenon, can be treated from eq. (26) 



If one defines the excitation times for a given trajectory by: 



i/*(X) = inf {X(t) e M} (27) 

t>«/*_i(*) 



with vq — 70 = 1, the energy configuration at time T, for an initial condition X is: 

m(T,X) 

X(T)=DF£.X+ ^x^-e^x, ( 28 ) 

where m(T, X) is the number of excitations on a time interval of length T for the initial condition X. The first 
term corresponds to the redistribution of the initial energy configuration while the second one corresponds to the 
redistribution of the energy quantum 6 = 1 dropped in the system at times i^-(X). Since the equilibrium average is 
assumed to be independent of the initial condition, the first term has to decay to zero as t — > oo, with a decay rate 
corresponding to the characteristic relaxation time to equilibrium. 

It is therefore important to well understand the (spectral) properties of the DF t X in the infinite time limit. Were 
S(X) be the Laplace operator, were the spectrum of DF l X be composed by Fourier modes, and the relaxation time 
to equilibrium would be the slowest mode. However, the mere presence of a singular term Z(X) certainly makes a big 
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difference. Since S depends on X one clearly has to study the decay rates averaged on a full (typical) trajectory or 
equivalently to compute the ensemble average. In this view, the law of the stochastic process {Z (X.(t))}^f^ (namely 
the density of active sites and correlations at all times) certainly plays a role. 

The numbers characterizing the decay (resp. expansion) rates of the norm of a small pertubation in the trajectory's 
tangent space of a point X under the action of the infinite product matrix DF l ~K, t — > oo are the Lyapunov exponents. 
They are real numbers, well defined under some moderate assumptions on DF~K (see Jig]) and are almost-surely 
independent of X. Furthermore they are also independent of the norm (in finite dimension). 

As shown in p"l[ ] and widely discussed in this paper all the Lyapunov exponents are different from zero, for finite 
L (weak hyperbolicity) . One remarkable consequence is that the asymptotic dynamics lies onto a fractal attractor 
and that the Lyapunov spectrum is closely related to the (local) fractal properties of the invariant set through the 
Kaplan- Yorke and the Ledrappier- Young formula jl8|,[ll). At this point a remark is necessary. Hyperbolicity is 
clearly a particular feature of Zhang's model and of similar models where the amount of redistributed energy from a 
critical site i depends on its energy X{. By opposition, in models like BTW's or Dhar's model the amount of transfered 
energy is a constant. As a direct consequence, in these models, the dynamics is simply a piecewise translation in the 
phase space, and the uniform measure is preserved |i] Hence, all lyapunov exponents are zero ]l9[j . There is therefore 
clearly a structural difference between the dynamics of Zhang's type models and sandpiles. This observation seems 
at first sight to ruin the hope to classify the Zhang's type models and the sanpiles in the same "universality class" . 
However, we show in this paper that hyperbolicity of Zhang's model is partially lost in the thermodynamic limit. 
Namely, some of the Lyapunov exponents go to zero as L — > oo, with a polynomial rate (exponent t\ in the last 
section) closely related to SOC critical exponents. It might therefore well be that these two class of model share the 
same SOC critical exponents in the thermodynamic limit, though their dynamics are still of different nature, even in 
this limit. 

Due to the skew product structure, the tangent map at any point X admits a natural splitting DF = 
(n u (DF , ir s (DF x)) where the one dimensional map jr"(DF*) is expansive. Indeed, the average expansion rate is 
given by : 

\ L (0)= lim hog{det(Tz u {DF\)))=Q L .log{N) (29) 

since det(TT u (DFi)) = N^=i Therefore, since ljl ^ 0, there is & positive Lyapunov exponential the dynamics. 

Note that this is due to the excitation rule, and that it reflects simply the "chaotic" properties of the Bernoulli shift. 

A more important issue concerns TT s {DF^y) = DFy^. The Oseledec theorem Jl6| asserts that under mild conditions 
on -DFx there exists a hierarchy of Lyapunov exponents Al(1) > ... \l(N), Lebesgue almost-surely constant if \±l is 
the SRB measure, and a hierarchy of nested subspaces (Oseledec splitting): 

R N = V 1 {±)dV 2 (X)dV n {±) 

depending on X, such that the norm of a perturbation v € Vj(X) \ V,+i(X) is given by : 

H-DF^.vH =C(X,i)e A ^-*.||v|| (30) 
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where lirut_,oo jlogC(K, t) = almost surely, namely Ai(i) is the exponential rate of variation of ||v||. Define 
M(X,i) = DF^.DFx and A = lim^oo M (X, t)& ,(the Oseledec multiplicative ergodic theorem insures that this 
limit exists almost-surely and is a constant) . Then the Lyapunov exponents are the logarithm of the eigenvalues of 
A. M(X,£) being symmetric it admits an orthogonal basis {vj(X,t)}- -, and eigenvalues /ij(X,i) such that Al(i) = 
lim^oo irlog(p,i(lL, t)). Furthermore, each Vj(X,t) converges exponentially to a vector Vi(X) in IR^, depending on 



X 20 1 . The Vj(X)'s constitutes therefore a basis for the Oseledec splitting. We call them the Oseledec modes for 
the trajectory of X. They can be numerically obtained from the QR decomposition used in the computation of the 



Lyapunov spectrum (see |17|]). It has been shown in |11| that the Al(i) are all negative for finite L, namely all vectors 
in R N are asymptotically contracted. 

From this discussion, one expects a close connection between the Lyapunov spectrum and the energy transport in 
the Zhang's model. In particular, the following formula can be proved [p"l| 

N 

V X L {i) = log{e).(l - ^l) — — — — ■ = p L log{e). < s > L (31) 

U <t>l 

It relates the Lyapunov spectrum, which characterizes local properties of the microscopic dynamics, to the avalanche 
statistical properties of the macroscopic system. Note that the exponent Xjj(i) gives the contraction rate in the 
direction Vj(X) versus the local time. One can also define the average contraction per avalanche, xl('), given from 
eq. © by: 

Xl{%) = — XlW = 32) 

1 - W£ p L 

Then, the sum of xi's, giving the average volume contraction per avalanche, is related to the average avalanche 
size by : 

N 

Y,XL(i)=log(e).<s> L (33) 
i=i 

Note that < s >l corresponds to the total energy transport within one avalanche and is believed to be related to 
the total response function [ fL5| . Our formula shows that it is also equal to the volume contraction in phase space 
produced on average by one avalanche. 



B. Oseledec modes. 



To each negative Lyapunov exponent Al(i), i = 1...N is associated a characteristic time ti,{i) = \\t{i)\ \ the 
time for of a perturbation in the Oseledec direction i to vanish. This defines therefore a hierarchy: 

t L (l)>t L (2) >...t L (N) 

Note that there is no contradiction with the expected critical behaviour in the thermodynamic limit, since as L — » oo 
there are an infinite number of characteristic time scales. 

From the physical point of view a perturbation can be viewed as a small modification of the initial energy landscape 
X. It can be localized (for example one site perturbed) or spread. The attenuation is due to two distincts effects : 
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• Propagation through the lattice. 

• Dissipation at the boundaries. 

Note that according to the Oseledec mode under consideration the contraction can be due (on average) to the 
effect of one avalanche (if tt{i) is small compared to the average avalanche size), or to the cumulative effect of many 
avalanches (if £/,(«) is large). The coefficient Xh{i) (eq. (|32])) gives the average contraction per avalanche for the i-th 
Oseledec mode. Therefore the number — Vr gives an estimate of the number of avalanches needed to have a reduction 

Xl(i) b 

of the initial perturbation of a ratio ^ for this mode. Therefore a crossover point can be estimated by : 

Xl(i c ) ~ 1 (34) 
We will call slow (resp. fast) modes the Oseledec modes corresponding to A^(z) << Aj,(i G ) (resp. A/,(i) >> A^(i c )). 

1. Bounds on the first negative Lyapunov exponent. 

We give a bound on the first Lyapunov exponent related to the energy dissipation at the boundaries. Call 
$° ut (i,X) = 1 — X^iLi^^Xy - Since the energy is locally conserved, $J"*(t, X) is the ratio of the initial energy 
Xj given by the site j to the boundary 8 A in t time steps. In other words, the energy coming from Xj and lost at 
time t is X).A,. The following holds: 

Proposition 1 The largest negative Lyapunov exponent, A^(l) admits the following bounds: 

> lim -log(l - min($° ut (t, X))) > A L (1) > lim -log(l - max X)) (35) 

t — >oo t 3 t — *°° t 3 

This is interpreted as follows. As t — > oo, <&°"*(i,X) — > 1, Vj, VX, since, eventually, all the initial energy coming 
from X has been lost at the boundaries. The limit lim^oo jlog(l — &° ut (t, X)) gives the exponential rate of conver- 
gence of §° ut (t, X) to 1. In other words, it gives the exponential decrease for the ratio of the initial energy still in the 
lattice at a given time. The maximal negative Lyapunov exponent is bounded by the extremal dissipation rates. One 
sees therefore that the contraction in the principal Oseledec mode is mainly due to the dissipation at the boundaries. 
We shall see later on that Al(1) is essentially related to the so-called dissipation rate. 

Proof It is easy to show that there exists a time t s depending on E c , e, d such that, whatever X each site in the lattice 
has relaxed at least once after this time and therefore all sites of the boundary have dissipated energy. At time t the 
energy coming from a site j with initial energy Xj and redistributed into the lattice is Yli=i -^-^x ij-^j- F° r t > t s , 
$° ut (t,X) > and J^iLi ij ^ s bounded away from 1. Since DF^ is a matrix with positive entries: 

JV N 

min(V GFt,,) = 1 - max $>? ut (t, X) < pfDF^) < max(V DF^ u ) = 1 - min m ut (t, X) < 1 (36) 

3 ~ f 3 3 r-f 3 

where p(DF t x ) is the spectral radius of DF^. 

The largest negative Lyapunov exponent is given by : 
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X L (1) = lim -log(\\DF^\\ 2 ) (37) 

t— *oo t 

where || H2 is the Li norm. In eq. ( |37| ) the limit does not depend on X, provided X belongs to the support of /t£. 
One has p(-DF^) < ||-DF^|| 2 and therefore: 

Al(1) > lim -logCl - max X)) 
' j 



Furthermore, all norms being equivalent in finite dimension eq. ( p7| ) holds also for the L\ norm where ||Z)i^||i = 
DF x.,ij \ X j\ 

Therefore, 



sup 1 • r ^^ ie DF Xi j's being positive the supremum is certainly achieved for positive Xj values. 



||DJ^l| 1 = sup ^ =U d =l-inf ^^(t.X).^ < l-infmin($f'(t,X)) 

The limit limr^oo log(l — min^ (<i>° nt (t, X))) is a constant for any X in the support of p.^. Hence: 

A L (1)< lim i/o. 9 (l-min($° M *(t,X))) 



□ 



2. Stabilizing modes 



The contraction in the principal Oseledec mode (first negative Lyapunov exponent) is mainly due to the dissipation 
at the boundaries. On the other hand, it is possible to have a large contraction in one local time step without 
reaching the boundaries. Indeed, the tangent matrix DF& has the following property, which can be checked by direct 
computation. 

Proposition 2 Let A = A C (X)©A„(X) where A C (X) = {i E A | Xi > X c }, A„(X) = {i € A | Xi < X c } andn c (X) = 
#A C (X), then DF^ has n c (X) eigenvalues e corresponding to the eigenvectors 

fci(x) = 2.d. ei - ^ e i ; ie A c( x ) (38) 

where Ui denotes the set of sites in A at distance 1 from i. There are also N — n c (X) neutral eigenvalues associated 
with the eigenvectors e», % S A n (X). 

The eigenvectors ki produces arbitrary large contraction as e — > 0. In particular, in the original Zhang's model 
where e = they correspond to kernel modes, which have eigenvalues 0. Note that, in this case, the dimension of 
the kernel of the product tangent map DF^ increases with t. However, it is strictly lower than N as t — > 00 pl|| . It 
is easy to check that these modes have zero energy, except if some of the e;'s correspond to sites neighbouring the 
boundary. This occurs with small (but non zero) probability. These modes act as directions where a single local time 
step is sufficient to reduce the initial perturbation by a factor e, with small variation of the total energy on average. 
They dynamically correspond to directions transverse to the attractor, and their contraction corresponds to a fast 
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convergence onto the attractor. For this reason we call them stabilizing modes. In the Lyapunov spectrum they can 
be identified because the corresponding Lyapunov exponents go to — oo as e — > while the other part of the spectrum 
weakly depends on e (see Fig [l]) . 




epsi!on=0.1 
epsiion=1 e-7 



FIG. 1. Lyapunov spectrum for e = 0.1 and e = 10 7 , L = 15, E c = 2.2. 



3. The transport operator. 



It is usually not possible to give an explicit formula for the whole Lyapunov spectrum, except in some specific cases 
PH . We propose here a mean-field ansatz which gives good results for the slowest modes, and has a nice interpretation 
in terms of random walk. It is based on the following observations. 



The Lyapunov exponents are the eigenvalues of the matrix A = El [A] 



Er 



linit 



(df^.df^ 



Since the matrix DF^.DF^ is bounded in L2 norm, VX,Vt, one has from the Lebesgue theorem : A 



lim t _ 00 E L 



On the other hand, the matrix 

C(t) = E L [DFx] 



(39) 



characterizes the (ensemble) average energy transport in t time steps. However, the connection between C(t) and A 
is loose. 

Were the transport be normal, namely were -DFx be independent of X and of the form £>Fx = I + 7-A, where 7 
is some constant, then were C(t) be equal to (/ + 7. A) . In this case, C(t) would be constant and symmetric. Then 



E L 



= I + 7. A and A = / + 7. A. In this case the Lyapunov exponents would be the eigenvalues of 
a one-step transport operator L = I + 7. A (Fourier modes). 

More generally, the (naive) hope would be that of finding an effective transport operator £ such that C(t) = £* 
and whose singular values (or eigenvalues if £ is self-adjoined) would give the Lyapunov spectrum. There is however 
a priori no hope for finding such an operator, in general. Note in particular that the assumption of independence 
of the matrices DF X (t) , first step towards a mean-field approach, is not a sufficient condition. Since, in this case 
El [UFjf] — El [f Fxf, one is lead to propose C — El [DFx] as a one-step operator. However, one needs further 



conditions to insure that the singular values of L give the Lyapunov exponents (see for example pi 22|). It appears 
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nevertheless that in the Zhang's model an effective transport operator can be found from a mean-field treatment 
which well approximates the slowest modes. 

The first obstacle towards a mean-field approach lies in the independence assumption. The matrix C(t) 
is a sum of time correlations terms of the form El [iS(X(t r _i))iS(X(i r _2)) • • ■ £(X(£o))] whose entry (i,j) 
wri tes Efa,...,*^ • • • A; 2)il A iuj Prob [Z ir _, (X(i r _0) = 1, . . . , Z h (X(*i)) - l,Zj(X(i )) = l]. Clearly, the 

non-vanishing terms in this sum are those corresponding to a path from j to i where each intermediate site has 
been active at least once with a non zero probability. A simple glance at this formula shows that a priori all time cor- 
relation functions of the joint probability of active sites, Prob \Zi r _ 1 (X(f r _i)) = 1, Zj 1 (X(ti)) = 1, Zj(X(io)) = l] 
have to be considered. 

However, Zhang's model, as a hyperbolic dynamical system, has exponential correlation decay (for finite L). The 
largest correlation decay rate is given by the entropy (VL-log(N). This decay rate is quite faster than the char- 
acteristic times related to the slow modes (for example the correlation decay rate of a site with itself is about 
—0.025 for E c = 2.2, e = 0.1, L — 20, corresponding roughly to the 320-th exponent in the spectrum, while the 
slowest Lyapunov exponent value is —0.000209871). More generally, we show in the last section, that Al(1) ~ jj 
to be compared to the decay rate u)L-log(N). On the other hand, for the slow modes, a small perturbation has 
essentially no variation during one step of an avalanche. In other words, the slowest Oseledec modes are not 
sensitive to the fast changes (one local time step) of the individual matrices DF^ (resp. the fluctuations of the 
Zj(X.(t))'s) but, rather, to the average variations on the characteristic time scale ii(«) = \~m > which is quite 
longer than a local time step. This suggests that one should consider the projection of the matrices ZJFxm's 
on the slow Oseledec space as independents. This leads to propose El [-DFx] = I + a.A.pL-I as a one-step 
transport operator. Note that we obtain the same result by assuming that the .Zi(X(i))'s are independent. In- 
deed, in this case Prob [Z ir _, (X(i r _i)) = 1, . . . , ^(Xfa)) = 1, ^-(X(i )) = l] = p L (ir-i) ■ ■ ■ Pl^PlU)- Then 
£(*) = ELi CH(&.PL.I) k = (I + aA.pL.iy. 

This approximation gives correct results .... provided that one multiplies the density of active sites by 2 !!! This 
approximation neglects indeed an important effect. Provided E c > a site cannot relax two successive time steps 
jn]], and therefore it relaxes at most only half of the time during one avalanche. This means in particular that the 
random variables Zi(~X.(t)), Zi(X.(t + 1)) are not independent and that the probability that one site relaxes at time 
t + 1 depends on what happened at time t. Besides, two neighbouring sites cannot be simultaneously active. In 
a certain sense, the lattice is "blinking" : during one avalanche all active sites are at pairwise distance |TT| . This 
therefore introduces strong correlations between Z(X(i)) and Z(X(f + 1)). 

One can, however, circumvent the problem by reparametrizing the time and considering the evolution of the 
process every two times steps. Equivalently, one replaces the stochastic process {Z(X(t))}i^ by a new process 
{Y(t')} t ^IJ L = {Z(X(i)), Z(X(t + whose components Yfc(i) take values in {0, l} 2 , where the event (1, 1) has 

zero probability and where t' = |. One can then encode the Yfe(i') values by — > (0,0) (no relaxation at times 
t,t + 1) and 1 — > (0, 1), (1, 0) (relaxation at time t or at time t + 1). This leads to the definition of a new "density of 
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active sites" p' L (i) = Prob [Y^t') = 1] = Prob [2,-(X(f)) = 1 or Z,-(X(i + 1)) = 1]. Since the events {^(X(t) = 1} and 
{Zi(X(< + 1) = 1} are disjoints, we have: p' L (i) = Prob [Z(X(i)) = 1] + Prob [Z(X(* + 1)) = 1] = 2.p L (i). Assuming 
now that the Yfc(i')'s are independent and considering p' L (i) as the effective density of active sites, one obtains an 
effective transport operator: 

£ = I + 2aA.p L .I (40) 
Calling 7i the singular values of £, our mean-field ansatz suggests that the slowest Lyapunov modes are given by : 

X L (i) = log( 7i ) (41) 

Note that this operator is self-adjoined for the metric pL-I and that the corresponding matrix can be made symmetric 

_ i 

by the variable change p L 2 .1. 

To check the validity of this ansatz, we first computed the density of active sites on a 20 x 20 lattice and found 
numerically the 7^'s from these data []. At the same time we computed the Lyapunov spectrum. A plot of the two 
curves is drawn fig. ^. One finds a very good agreement over a large part of the spectrum, and the discrepancy 
increases towards small times scale, as expected. 
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We weren't able to go beyond L = 20 in the Lyapunov spectrum computation. We used a version of the Eckmann-Ruelle 
algorithm Jl7| revisited from Von Bremmen et al. p3[ |. Nevertheless, we needed two weeks of computation on a Pentium II 300 
for the case L = 20, with a relative accuracy of 10 . 
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FIG. 2. Lyapunov spectrum, logarithm of the C singular values, and normal diffusion modes for E c — 2.2,e = 0.1, L — 20. 
Fig ^ : 80 first modes; Fig : Full spectra drawn with lines in order to better see the shape. 

4- The role of the spatial variations in the density of active sites. 

It is usually assumed in the SOC literature, when dealing with the model's spatial properties, that only the density of 
active sites, and more precisely, its lattice average p^f, has to be taken into account. One therefore neglects the spatial 
dependence of pr,. In our approach this would lead to an effective transport operator / + 2.apf^ A, corresponding to 
normal transport. In this case, the slowest Oseledec modes would simply be the eigenmodes of the Laplacian with 
zero boundaries conditions on 9 A, and the Lyapunov exponents would then correspond to the normal diffusion modes: 

d 

\ L (i) = log{\l + 2.a.pf '•E( C0S (St ) ~ 1)D (42) 

fc=i 

where the Laplacian modes are parametrized by the quantum numbers n = (ri)-),k = 1 . . . d, sorted such that the 
corresponding eigenvalues are decreasing (and i refers to the placement of the exponent in this sequence). In fig. 

, we also plot the diffusion eigenvalues of eq. (fl2|). The computed Lyapunov exponents are different from these 
values except for the largest ones. This approximation is therefore too crude and gives a wrong spectrum. Note in 
particular that the shape of the spectrum differs, namely the discrepancy cannot be corrected by a mere multiplication 
of p"£ by some factor. Since the Lyapunov exponents contain all the relevant informations about the dynamics at 
stationarity, our conclusion is that the non- homogeneity of plays a key role in computing dynamical quantities, 
and implies unfortunately that the zero-th order "mean- field" approaches, which approximate the density of active 
sites as a constant lead to incorrect estimates for finite size when dealing with intermediate time scales. On the other 
hand, this should lead to correct results when dealing with the longest time scales, since the first modes are well 
approximated by a transport operator where ph(i) is considered as uniform. 

In the litterature one often meets an (apparent) contradiction (see for example the original paper form Zhang J5j and 
subsequent analysis by Pietronero et al. p4j ) . One assumes the transfer of energy on large time scales to be normal 
diffusion, while in the same time an anomalous diffusion exponent z ^ 2 is computed. It was certainly clear in the 
spirit of these authors that one has to distinguish the average transport on many avalanches (long time scales) from 
the average transport within one avalanche (characterized by z). Our results on the Lyapunov spectrum makes this 
distinction quantitative. We show that the transport on the longest time scales is normal with a good approximation, 
while the average transport within one avalanche, namely on time scales corresponding roughly to the crossover point 
Az,(i c ), is clearly anomalous. 

5. The random walk picture. 

The operator C is the Laplace-Beltrami operator associated to a random diffusion in a "medium" or a landscape that 
is not flat, corresponding to the metric g — pl-Ij where / is the identity matrix on IR W . It has a nice interpretation in 
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the so-called random walk picture []. Assume for a moment that the energy of a site is composed by (undivisible) energy 
quanta r\ that can be made arbitrarily small (this is way to "map" the Zhang's model to a sandpile) . Assume that we 
are in the stationary regim, and that at the initial time we drop a grain in some place and study its motion. At each 
time step where it is involved in a relaxation process this grain makes a jump at random in one of the 2.d directions 
in the lattice. From this point of view, the stochastic dynamics of the grain is driven by the underlying dynamics of 
eq.(Q). If we assume that the evolution is Markovian, the probability of jumping from i to some nearest neighbour 
j depends only on the state of i at time t and is given by a transition rate Wij — a.pi,{i), while the probability of 
staying at the same place is 1 — pi,(i)(l — e) (remember that only an amount a = of the energy is transferred to 
the neighbours when a site relaxes). From this consideration, one obtains the equation for the probability P(i,t) that 
a grain is at place i at time t, before it leaves the lattice P(., t + 1) = [I + a.Ap£,].P(., t) = [I + q;.A/5l]*.P(., 1), and 
one recovers the operator obtained above when assuming that the Zj(£)'s were independent. Indeed, the independence 
assumption of the Zj(i)'s is equivalent to the Markovian assumption in the random walk picture. The probability of 
jumping for a grain at time t depends a priori on its whole past via a Chapman-Kolmogorov equation whose transfert 
matrix is a sum of terms containing conditional probabilities 

ProblZ^iXiU^)) = l|Z 4r _ 2 (X(i r _ 2 )) = 1, . . . , Z^Xih)) = l,^.(X(t )) = l] = 

Prob [^(Xfc-i)) = l,Z ir _ a (X(t r - 2 )) = 1, . . . ,Z n (X(Q) = l,Z 3 -(X(t )) = l] = 
Prob[Z ir _ 2 (X(t r - 2 )) = l,...,Z 4l (X(t!)) - l,^(X(t )) = 1] 

Prob [^(Xftr-i)) = 1] 

where the last equality holds when the Zk(t)'s are independent. In this case, Wij = a. Prob [Zi(X.(t)) = 1] = api,(i) 
for the j nearest neighbours of i. 

However, we saw above that the process is not strictly Markovian since a jump from a given site cannot take place 
at two successive time steps. In other words, the probability of a jump i — > j depends on the state of i at time t 
and at time t — 1. The system has some memory (at least two-timc-stcps) . However, defining the random variables 
Yfe(t)'s, as above, and assuming them to be independent amounts to render the random walk Markovian by a suitable 
reparametrization of the process, and gives a transfert equation 

P(., t + 1) = [/ + 2.a.Ap L ].P(., t) = CKP{., 1) (43) 

Therefore, the operator C characterizing the decay of a small perturbation can also be interpreted as the transfert 
matrix of a random walk in a medium where the diffusion rate depends on the location. 



7 B.C. is very grateful to P. Grassberger and D. Dhar for illuminating discussions on this point in Trieste. 
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6. Density of active sites and average energy flowing towards the boundaries. 



Eq. ( p3| ) characterizes the energy transport in the lattice, but does not take into account the source term (addition 
of a grain) required to reach stationarity. Indeed, each times a grain exits the lattice, one must add another grain at a 
random place i, with probability U)i,(i) ; this is a source term. Call Vl the equilibrium state of the random walk and 
V(d\) the set of sites at distance one from the boundary. The probability for a grain to exit is 2a X^ev(dA) PL^t-^hii)- 
It is obviously proportional to the outgoing energy flux, which is, at stationarity, equal to the incoming flux (resp. 
the probability of adding a grain in the lattice) , namely 0: 

Q L = 2a ]T p L (j)-X L (j) (44) 

The complete equation for the energy at stationarity is : 

2a.A[p L .X L }+Cu L {i) =0 (45) 

with zero boundaries conditions and with the constraint (|44|). 

In this equation one distinguishes a local transport term, and a source term which depends on a global constraint. 
When the excitation is uniform eq. ([l5|) reduces to A[p^.Xi] + 2 ^ L L d = 0. 

The difficulty in solving this equation is that it deals with the product p^.X^. On the other hand, it is known in 



the literature that X/, converges to a uniform energy distribution over the lattice as L — > oo 25 1. Assume now that 
we can write X^ as: 

X L = X + f(L) (46) 

where ||/(L)|| goes to zero as L — > oo and where Xo is spatially uniform, i.e. Xo(i) = const = xl- At the zero-th 
order, one obtains for pr, the following equation: 

where L d XL = E tot , the average total energy in the lattice. The solution of this equation can be easily found by 
decomposition on the cigenmodes of the Laplacian. The general solution is: 

d 

p L (x) = ^2A a J[sin(ki.Xi) (48) 

n i— 1 

where n = (n±, . . . , rid) is the set of quantum numbers parametrizing the eigenmodes of the discrete Laplace operator, 
s n = 2(y] i -_ 1 cos(ki) — d) is the corresponding eigenvalue with ki — jxj, 

A 2*-\u> L nti^ 



aE tot (L + l)° 



The cautious reader has noticed that this equation is not dimensionally correct, since no energy term appears on the l.h.s.. 
One should indeed multiply by S, the input energy quantum, which is set to one throughout this paper. 
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and 

x=l Sln \2(L + l)) 

where = 2m^ + 1. Surprisingly, this already gives quite a good approximation for p^, which becomes better and 
better as L increases (see fig. and fig. in the next section.). 




Away from the boundaries, one expects rotational invariance for p^(x). This can be checked by expanding the 



function sin near to Xi 



1 ... d up to third order. One obtains the well known paraboloid form p6j pl(x) 



if — Ki J2i=i x \i where the constants K , K\ can be easily deduced from eq. ( fi7| ) 
One also obtains the average density of active sites, p°£ — -he J2iLi 



av 

Pl 



E 



(49) 



L d (L + l) d aE tot 

which is expected to hold for sufficiently large L. We give a plot fig. [| where it clearly appears that this formula gives 
already a quite good estimate for L = 15. 



Density of critical sites lattice average 
First order approximation 



FIG. 4. Plot of pT and solution of eq. @) versus L, for E c = 2.2,e = 0.1. 
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IV. SCALING PROPERTIES OF THE LYAPUNOV SPECTRUM. 



Zhang's model, as a hyperbolic dynamical system, cannot exhibit a critical behaviour for finite size, since it has 
exponential correlation decay []. However, since a critical behaviour is conjectured in the thermodynamic limit, one 
expects that hyperbolicity is lost as L — > oo, namely some of the Lyapunov exponents go to zero. It is therefore of 
crucial importance to know the behaviour of the Lyapunov exponents as L — > oo. In this section, we first discuss the 
time scale separation between activation rate and dissipation rate, which is believed to be a fundamental ingredient 
to have SOC, and its links to the Lyapunov exponents. We then show that using a Finite Size Scaling ansatz provides 
a scaling exponent from which the scaling of some SOC observables can be obtained. 

A. Time scale separation. 



L a ah t ot 



Equation (|49f) can be written as: 



TT d c 2 

where (L + l) d -jL — J2 n s * — ~" ^et us estimate the scaling of this sum as L — > oo. First, set d = 1 and fix a > 
arbitrarily small. The sum over n — n\ can be split into a part such that n < (L + l)a and another part where 
n > (L + l)a. In the first sum, C n ~ 2 ^^ 1 - ) , s n ~ — ^ J +1) 2 ' wm le the second sum is smaller than (L + l)C(a), where 
C(a) is bounded for a > 0. Therefore jl ~ {L + l) 3 .^, where S ~ Y^ n <(L+i)a n* stays bounded as L — ► oo. Then 
7l ~ (L + l) 3 . This argument can be generalized for any d by splitting the sum over n = (m, ri2, ■ ■ • Ud) into sums 
where k indexes are smaller that a(L + 1), k going from to d. It is easy to see that the main contribution is due 
to the terms such that d indexes are < a(L + 1), giving a leading contribution ^2(L + \~) 2d+2 and 7l ~ L d+2 . We 
therefore conclude that p°£ scales like : 

q l (L + i) d + 2 ^J? = ^J? 

lL E tot L d E tot L d X 1 ' 

Set h = yjt for the driving rate and call e = -fer = xl.L~ 2 . One obtains the energy conservation equation : 
h = p% v e and therefore e is the energy dissipated per active site and per unit time. It corresponds to the dissipation 
rate introduced by Vespignani et al. flTq] . Since < xl < E C ,VL, xl plays no role in the asymptotic scalings in L 
and therefore e ~ L~ 2 , as already anticipated by a mean-field approach in [T^ ]. 

The average value of observables like size, duration, etc ... is known to diverge with a power law scaling (x) L ~ L lx . 
Therefore j~y- — > like L~ 1t as L — > oo where y T > 1 Since < iDl < 1 (see (pC|)), eq. ( pl| ) implies that 
Pl = f(jpj-) = — y^r + O(^g-). This is in particular clear for i? c < 1, since = iD^, which implies 
Pl = i+{ T ) an< i therefore ai = 1, a 2 = 1. For general _E C using this form gives, from eq. a! = 1 and : 



1 «2 



TTT" " TIT (52) 



9 The exponential correlation decay is a general property of hyperbolic systems but in presence of singularities one can also 
observe a polynomial correlation decay and weak initial condition senstivity M. 
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and 



as L — > oo. It follows therefore that 



We have therefore shown that 



° L ~ W L (53) 



Pf ~ L -~(r+2-d ^ 



h -> 0, e -> 0, p2 u = ~ -> 0, as L ^ oo (55) 

In |l5| Vespignani et al. discussed the necessity of this triple limit in order to have SOC. In their analysis the 
activation and dissipation rate where however free parameters (tunable "by hand" ) . In Zhang's model, h and e are 
not free, since they are fully determined by the dynamics. Therefore, we have shown that the three limits discussed 
by Vespignani et al. |J] are indeed achieved, without external fine tuning of some parameter, in Zhang's model, by 
the simple constraints one imposes on the dynamics (adiabatic driving). 

From ( |53| ) we have that the positive Lyapunov exponent (resp. the entropy) u)L-log(N) — > in the thermodynamic 
limit. On the other hand, the first negative Lyapunov exponent is given with a good accuracy by the normal diffusion 
operator 1 + 2.p a L v A (see fig.|), which implies that X L (l) ~ pTL^ 2 . Another way of arguing is to note that from 
theorem 1, Aj,(l) scales like the average ratio of energy dissipated by one site. From the local conservation of 
energy, Xl(1).Xl ~ h, then Al(1) ~ p1"L~ 2 . Therefore, Al(1) — > in the thermodynamic limit. Actually, the 
Finite-Size Scaling analysis of the next section suggests that a large number of negative Lyapunov exponents also 
go to zero as L — > oo. But the double limit A^(0) — > 0, A^(l) — > already show that the hyperbolicity is lost in 
the thermodynamic limit. Note however that these two exponents are not indepedent since local conservation of 
energy imposes Al(1).5l <~ h which implies ~ log(N) • Incidentally, this validates the separation of time scale 

between the correlation decay time j^qj and the largest transport characteristic time T~rn we use d when deriving 
the mean-field transport equation for the slowest modes, in the previous section. 



B. Finite Size-Scaling of the Lyapunov spectrum. 



An approximate expression for the modes related to the transport in the lattice is obtained from the operator C 
(eq. ^), whereas an approximate equation for p^ is given by eq. ([f9|). However, at the moment we don't have an 
analytic expression for the modes of C. In this sequel, we restrict to the scaling of the slowest singular values of C 
with the system size. 

When dealing with scaling analysis in the thermodynamic limit, one usually first tries to use Finite Size Scaling 
(FSS). This is a standard tool in statistical mechanics. It has also be proposed in SOC as an anstaz for the scaling 
of the probability distribution of avalanches observables 0. However, its validity has recently been asked in this case 
§• 
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Nevertheless, since this is certainly the first anstaz one can try to do scaling analysis, we try in this section a 
Finite-Size Scaling ansatz for the Lyapunov spectrum and look at the results and conclusions we are lead to [5 We 
assume therefore that, for any L, there exists a change of coordinates i — ► </>l(*),-^l ">Pl(Xl), depending on L, 
such that the points of the spectrum {i, Ai(-i)} are mapped onto the same "universal" curve ^ {x, A(x)}, where 
A(x) = ° o (f>~£ (x). Furthermore we assume (as in usual Finite Size Scaling) that the coordinate changes are 
simple dilatations where <Pl(x) = L@ x .x, iPl(x) = L /3x - Tx (x). Then: 



\{x) = L^- T \X L (x.L-^) 
Equivalently, knowing the curve {x, A(x)} the spectrum for a given size is 

X L (i) = L- ^X(iL-^), i = l...L d 
Since the set of indices i € {l . . . L d } it is evident that : 

Px = d 



(56) 



(57) 



(58) 



The exponent t\ can be numerically computed by several means. A first one is to minimize the euclidean distance 
between the spectra obtained for different lattice sizes, with respect to t\. Another way is to compute the sum of the 
Lyapunov exponents. Indeed: 



N 1 r l 

»=1 y=L- d L d 



X(y).dy 



Assuming that A(y) is bounded as y — * and that < K = f Q X{y)dy < oo one obtains: 

S L - K.L d - {1 - T ^ 



(59) 



which allows one to compute T\. The value of t\ for d — 2, e = 0.1 and different E c values are given Table 1. These 
values were obtained for a sample of spectra from L = 10 to L = 20. We note in particular that t\ depends on 
E c . At the moment we have no way of knowing wether this effects persists in the thermodynamic limit. Note that 
these values are given as indications but the correct estimation of t\ certainly requires further investigations for 
consequently larger system size. These numerical studies are beyond our present computer performances. 



Ec 




0.6 


0.632 


1.1 


0.622 


1.5 


0.621 


2.2 


0.560 


4.1 


0.524 



Table 1: Computed values of t\ versus E c , obtained from eq. (^), for samples of size L=10 to L=20. 



Note that FSS of the Lyapunov spectrum is not a general property of dynamical systems, even close to a phase transition 

point nn 

11 Note that this curve depends on the parameters E c , e, d 



23 



The data collapse of spectra is drawn Fig. |5|. Though a good data collapse is not sufficient to ensure FSS, Fig. 
H indicates that this gives a good approximation of the spectrum. Actually, we don't expect the FSS to hold for 
the whole spectrum (in particular the kernel modes could have a different scaling). For the following discussion it 
is however sufficient to assume that FSS holds for the slowest modes. This is a reasonnable assumption since these 
modes are well approximated by normal diffusion. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



FIG. 5. Data collapse of the Lyapunov spectrum for Ec=2.2,e = 0.1,L=12,14,16,18. 

We now relate the exponents 7 S ,7 T (and other characteristic exponents like z, the anomalous diffusion exponents) 
to t\. Note that ^y x is related to the critical exponent t x and therefore our discussion suggests that there is a link 
between t\ and the critical exponents t s ,t t . 

The FSS leads to Al(1) = L~ d Tx .\(L~ d ). However the analyticity properties of A near zero are not known. Assume 

that X(x) ~ x a ,x ~ where a may depend on d (seemingly a = 1 for d = 2). Then Xl(1) ~ £~ d - TA ~ d ". From 
Al(1) ~ p1"-L~ 2 one obtains: 

p a L v ~ L- ±Tx+2 - da (60) 
and from (|4j) 7 r + d — d(r x + a). Finally, from eq. (||l]) , (|2|) , (||) one gets: 

dr\ = d - 7 S + 7 T (61) 

which gives: 

7, = 2 (62) 

~f T =dT X + 2-d (63) 

12 Under the finite-size scaling assumption of Pl(x) one finds that j x = (3 X .(2 — r x ) where L 13 " is the scaiing for the maximal 
value of a; in a lattice of size L. 
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and 

a = - d (64) 

The equation for 7^ has been already anticipated by many authors on the basis of numerical simulations [p5| , 
mean- field approach Jl5| and has been proved in Dhar's model for d = 2 by Dhar himself p9| . The equation for a is 
well verified in d = 1 and d = 2. This relation deserves however further investigation in larger dimensions. It suggests 
in particular that the curve \(x) is not C 1 at zero for d > 2, i.e. the largest exponents do not go to zero in a smooth 
way as L — > 00. 

Finally, the anomalous diffusion exponent z, characterizing the average transport within one avalanche, is equal to 
7r if one assumes that the average avalanche radius scales like L in any dimension Q|. Equivalently, one can notice 
that the crossover point for the xl(*)' s spectrum (eq. (|32|)) is ~ < ^ > and does not depend on L. From eq. ( |63| ) it 
follows that the transport on time scales of order < r >l is anomalous (z < 2) iff t\ < 1. Note however that this 
argument assumes that the FSS is still valid at the crossover point. 

This result suggests therefore that some of the critical exponents of SOC can be obtained from simple scaling ansatz 
on the Lyapunov spectrum. 



As a final remark, note that the E c dependance appearing in table 1 would have to be clarified since it suggests that 
the critical exponents depend on E c . This was already argued in Js| — fl3j] and suggested from numerical simulations 
(though not discussed) in pc| . Note, however, that the dependence of dynamical quantities in the control parameter 
in a dynamical system is more a rule that an exception. One certainly needs very special properties to ensure that 
the critical exponents are constant in the limit L — > 00, whatever E c . If this happened to be true it would mean that 
Zhang's model is somehow non-generic, at least from the dynamical systems point of view. 



V. CONCLUSION 



In this paper, we investigated the dynamics of Zhang's model in terms of the Lyapunov exponents and Oseledec 
modes. Due to the piecewise affine structure of the model, the Lyapunov exponents, usually related to the local 
properties of the dynamics (expansion rates, fractal dimensions, entropy), also appear as characteristic rates of energy 
transport in the system. We showed that the spectrum is roughly divided into two parts, the slow modes corresponding 
to transport and dissipation and the fast ones essentially associated with the stability of the attractor. Even if the 
Oseledec modes are the analogous of Fourier modes in normal diffusion, they are not normal, because the density of 
active sites is not spatially homogeneous. The slow Oseledec modes correspond rather to a diffusion in a metric which 
is not flat and given by the density of active sites. Only for the slowest mode are the Lyapunov exponents the same 
as for the largest rate in normal diffusion. This is important since the slowest mode characterizes the equilibrium 
properties of the model. This means that the usual mean-field approaches, which replace the density of active sites 
by its lattice average, are correct if one considers properties related to the longest time scales. Since the critical 
exponents 7 S ,7 T characterize statistical properties on the largest time scale, they are naturally related to the slowest 
Lyapunov exponent. 
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We investigated the scaling properties of the spectrum with respect to the lattice size and found that Finite Size 
Scaling gives a good approximation. In particular we extracted a critical exponent t\ which is related to the usual 
critical exponents computed in the litterature. However, there is clearly a lot more information in the Lyapunov 
spectra than in the usual critical exponents. 

The scaling form shows also that in the thermodynamic limit a part of the spectrum goes to zero, corresponding to 
translation invariance, and zero dissipation. In this way the Zhang's model is not hyperbolic in the thermodynamic 
limit. This limit has now to be studied in more details, especially as far as the vanishing of correlations is concerned. 
It may indeed be a way to make a connection between SOC and the usual critical phenomena. 
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